/*  primetbl.cpp  by K.Tsuru  */
// since ver 2.181
#ifndef TOOLS_H
#include "tools.h"
#endif  // TOOLS_H
/***************************
The number of prime under N.
****************************/
#if 0
ulong PrimeNumberUnder(ulong N){ // slow N > 10^8
    bool isPrime = false;
    long i, j, count = 2; // 2 and 3
    // start at 5
    for (i = 5; i <= N; i += 2) {
        for (j = 2; j * j <= i; j++) {
            if (i % j == 0) {
                isPrime = false;
                break;
            }
            isPrime = true;
        }
        if (isPrime) count++;
    }
    return count;
}
#else
/********************************************************
The rough estimate value for the number of prime under N.
pi(x) ~ x/ln(x) (Prime number theorem)
*********************************************************/
ulong PrimeNumberUnder(ulong N){
    double f, L = log10((double)N);
    int m = L+0.01;// to integer
    if(m <= 5) f=1.17;
    else if(m <= 6) f=1.11;
    else if(m <= 8) f = 1.073;
    else if(m <=10) f = 1.055;
    else if(m <=12) f = 1.045;
    else f= 1.04;
    double p = log(N);
    double M = f*double(N)/p;
    return (ulong)M;
}
#endif
/****************************************
It makes the prime table between 2 and n
using the sieve of Eratosthenes.
*****************************************/
ulong MakePrimeTable(NCBlock <ulong>& table, const ulong n) {
    ulong count = 0;
    ulong M = (n - 3)/2;    // from 2M+3 <= n see (*) below
    while(2*M + 3 > n) M--;
    bool* flag = new bool[M+1]; // sizeof(bool) = 1 (byte)
    ulong ts = PrimeNumberUnder(n); // table size

    table.size(ts+1, 0);
    table[count] = 2;   count++;
    ulong i, k, p;

    for (i = 0; i <= M; i++) flag[i] = true;

    for (i = 0; i <= M; i++) {
        if (flag[i]) {
            p = 2 * i + 3; // (*) p<n
            table[count] = p;
            for (k = i + p; k <= M; k += p) flag[k] = false;
            count++;
        }
    }
    table[count] = 0;   // set last element = 0
    delete[] flag;
    return count;
}